% Replication of the trivariate VAR in Stock and Watson (2001, JEP).
% Figure 1 and Table 1.B.
%
% The VAR Toolbox 2.0 is required to run this code. To get the 
% latest version of the toolboxes visit: 
% 
%       https://sites.google.com/site/ambropo/MatlabCodes
% 
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com


%% PRELIMINARIES
% =======================================================================
%clear all; clear session; close all; clc
%warning off all



%% PRELIMINARIES
% =======================================================================
%clear all; clear session; close all; clc
%warning off all

%_TFP specification_* %
tfp_sp=0;
% 0: no TFP data
% else(1): with TFP




% Load data
%[xlsdata, xlstext] = xlsread('SW2001_Data.xlsx','Sheet1');
%[xlsdata, xlstext] = xlsread('VAR_data_tes1.xlsx','hours_tfp');
%[xlsdata, xlstext] = xlsread('VAR_data.xlsx','hours_noA');
% Define and transform (if needed)

load intvar_monte_data;
macro=intvar_monte_data;
X = macro;
%vnames = xlstext(1,2:end);


%% VAR ESTIMATION
% =======================================================================
% Define number of variables and of observations
[nobs, nvar] = size(X);
% Set deterministics for the VAR
det = 1;
% Set number of nlags
nlags = 4;
% Estimate VAR
[VAR, VARopt] = VARmodel(X,nlags,det);
% Print estimation on screen
%VAR.vnames = vnames;
%VARprint(VAR,VARopt);
%disp(' ')
%junk = input('Press enter to continue');


%% COMPUTE IRF AND FEVD
% =======================================================================
% Set options some options for IRF calculation
VARopt.nsteps = 63;
VARopt.ident = 'oir';
%VARopt.vnames = vnames;
% Compute IRF
[IRF, VAR] = VARir(VAR,VARopt);
% Compute error bands
[IRFINF,IRFSUP,IRFMED, IRFBAR] = VARirband(VAR,VARopt);
%[pIRFINF,pIRFSUP,pIRFMED, pIRFBAR] = VARirband_price_level(VAR,VARopt);
% Plot
%VARirplot(IRFMED,VARopt,IRFINF,IRFSUP);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Creating IRF variables WITHOUT TFP%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if tfp_sp == 0
% 
% diffIRFINF_pai= [0;diff(IRFINF(1:62,3,5))];
% diffIRFSUP_pai= [0;diff(IRFSUP(1:62,3,5))];
% diffIRFBAR_pai= [0;diff(IRFBAR(1:62,3,5))];
% pIRFINF_pai=[0;pIRFINF(:,3,5)];
% pIRFSUP_pai=[0;pIRFSUP(:,3,5)];
% pIRFBAR_pai=[0;pIRFBAR(:,3,5)];
% 
% inf_mn=pIRFBAR_pai;
% inf_hb=pIRFSUP_pai;
% inf_lb=pIRFINF_pai;
% y_mn=IRFBAR(:,1,5);
% y_hb=IRFSUP(:,1,5);
% y_lb=IRFINF(:,1,5);
% i_mn=IRFBAR(:,5,5);
% i_hb=IRFSUP(:,5,5);
% i_lb=IRFINF(:,5,5);
% c_mn=IRFBAR(:,2,5);
% c_hb=IRFSUP(:,2,5);
% c_lb=IRFINF(:,2,5);
% %w_mn=IRFBAR(:,4,5);
% %w_hb=IRFSUP(:,4,5);
% %w_lb=IRFINF(:,4,5);
% h_mn=IRFBAR(:,4,5);
% h_hb=IRFSUP(:,4,5);
% h_lb=IRFINF(:,4,5);
% 
% 
% else
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Creating IRF variables WITH  TFP%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% diffIRFINF_pai= [0;diff(IRFINF(1:62,3,6))];
% diffIRFSUP_pai= [0;diff(IRFSUP(1:62,3,6))];
% diffIRFBAR_pai= [0;diff(IRFBAR(1:62,3,6))];
% pIRFINF_pai=[0;pIRFINF(:,3,6)];
% pIRFSUP_pai=[0;pIRFSUP(:,3,6)];
% pIRFBAR_pai=[0;pIRFBAR(:,3,6)];
% 
% 
% inf_mn=pIRFBAR_pai;
% inf_hb=pIRFSUP_pai;
% inf_lb=pIRFINF_pai;
% y_mn=IRFBAR(:,1,6);
% y_hb=IRFSUP(:,1,6);
% y_lb=IRFINF(:,1,6);
% i_mn=IRFBAR(:,6,6);
% i_hb=IRFSUP(:,6,6);
% i_lb=IRFINF(:,6,6);
% c_mn=IRFBAR(:,2,6);
% c_hb=IRFSUP(:,2,6);
% c_lb=IRFINF(:,2,6);
% %w_mn=IRFBAR(:,4,6);
% %w_hb=IRFSUP(:,4,6);
% %w_lb=IRFINF(:,4,6);
% h_mn=IRFBAR(:,4,6);
% h_hb=IRFSUP(:,4,6);
% h_lb=IRFINF(:,4,6);
% tfp_mn=IRFBAR(:,5,6);
% tfp_hb=IRFSUP(:,5,6);
% tfp_lb=IRFINF(:,5,6);
% 
% end





%figure(101)
%subplot(2,2,1)
%plot(1:24,IRFINF(:,3,7),':'); hold on;  plot(1:24,IRFSUP(:,3,7),':'); hold on; plot(1:24,IRFMED(:,3,7),'-') ; 
%title('INF');
%hold on

%subplot(2,2,2)
%plot(1:24,IRFINF(:,1,7),':'); hold on;  plot(1:24,IRFSUP(:,1,7),':'); hold on; plot(1:24,IRFMED(:,1,7),'-') ; 
%title('Y');
%hold on

%subplot(2,2,3)
%plot(1:24,IRFINF(:,2,7),':'); hold on;  plot(1:24,IRFSUP(:,2,7),':'); hold on; plot(1:24,IRFMED(:,2,7),'-') ; 
%title('C');
%hold on

%subplot(2,2,4)
%plot(1:24,IRFINF(:,7,7),':'); hold on;  plot(1:24,IRFSUP(:,7,7),':'); hold on; plot(1:24,IRFMED(:,7,7),'-') ; 
%title('i');
%hold on


% Compute FEVD
%[FEVD, VAR] = VARfevd(VAR,VARopt);
% Compute error bands
%[FEVDINF,FEVDSUP,FEVDMED] = VARfevdband(VAR,VARopt);
% Plot
%VARfevdplot(FEVDMED,VARopt,FEVDINF,FEVDSUP);


%% Print Table 1.B on screen
% =======================================================================
% Retrieve Forecast Error Variance Decomposition
%FEVD_Table(1, :) = FEVD(1,:,1);
%FEVD_Table(2, :) = FEVD(4,:,1);
%FEVD_Table(3, :) = FEVD(8,:,1);
%FEVD_Table(4, :) = FEVD(12,:,1);

%FEVD_Table(5, :) = FEVD(1,:,2);
%FEVD_Table(6, :) = FEVD(4,:,2);
%FEVD_Table(7, :) = FEVD(8,:,2);
%FEVD_Table(8, :) = FEVD(12,:,2);

%FEVD_Table(9, :) = FEVD(1,:,3);
%FEVD_Table(10,:) = FEVD(4,:,3);
%FEVD_Table(11,:) = FEVD(8,:,3);
%FEVD_Table(12,:) = FEVD(12,:,3);

% Print on screen
%disp(' ')
%disp('Variance Decomposition of Inflation (t=1,4,8,12)')
%disp('---------------------------------------------------')
%mprint(FEVD_Table(1:4,:));
%disp('Variance Decomposition of Unemployment (t=1,4,8,12)')
%disp('---------------------------------------------------')
%mprint(FEVD_Table(5:8,:));
%disp('Variance Decomposition of Fed Funds (t=1,4,8,12)')
%disp('---------------------------------------------------')
%mprint(FEVD_Table(9:12,:));